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SECTION  I 
INTRODUCTION 


Over  the  past  several  years,  there  has  been  considerable  effort  devoted  to  the 
development  of  mathematical  models  to  represent  the  nonlinear,  inelastic  behav- 
ior of  geologic  materials.  These  models  are  necessary  for  the  accurate  analy- 
sis and  prediction  of  ground  motion  response  and  structure/medium  interaction 
under  explosive,  earthquake,  and  vibration  excitation.  Several  material  models 
are  currently  in  use;  e.g.,  variable  moduli  models  and  various  elastic/plastic 
models  (ref.  1). 

The  objectives  of  this  effort  were  to  investigate  a particular  capped,  elastic/ 
plastic  material  model  and  to  implement  this  model  into  a wave-propagation  code 
for  use  in  ground  motion  calculations.  Weidlinger  Associates  have  made  signi- 
ficant advancements  in  the  development  of  cap  models  (re.''s.  2,  3).  The  work 
reported  here,  much  of  which  is  based  on  their  research,  was  undertaken  to  com- 
plement earlier  work  (ref.  2)  by  providing  a fairly  detailed  analysis  of  one 
specific  form  of  the  general  soil  cap  model. 

The  general  soil  cap  model  is  an  elastic,  nonideally  plastic  model  with  yield 
criteria  that  combine  both  ideal  plasticity  and  strain  hardening.  This  model 
has  been  used  to  simulate  the  uniaxial  strain,  standard  triaxial  stress,  and 
proportional  loading  tests  performed  in  laboratory  experiments  (ref.  1)  and  al- 
so satisfies  theoretical  requirements  for  uniqueness  and  stability.  A brief 
summary  of  the  general  soil  cap  model  and  a similar  model  for  rock  is  contained 
in  reference  1. 


1.  Nelson,  I.,  Baron,  M.  L.,  and  Sandler,  I.,  "Mathematical  Models  for  Geologi- 
cal Materials  for  Wave-Propagation  Studies,"  Shock  Waves  and  the  Mechanical 
Properties  of  Solids^  Syracuse  University  Press,  Syracuse,  New  York,  1971. 

2.  Sandler,  I.,  and  Rubin,  D.,  A Modular  Subroutine  for  the  Cap  Models  Report 
DNA-3875F,  Defense  Nuclear  Agency,  Washington,  D.C.,  January  1976. 

3.  Sandler,  I.  S.,  DiMaggio,  F.  L.,  and  Baladi,  G.  Y.,  "Generalized  Cap  Model 
for  Geologic  Materials,"  Journal  of  the  Geotechnical  Engineering  Division, 
ASCE,  Vol.  102,  No.  GT7,  July  1976,  pp.  683-699. 
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Subsequent  sections  of  this  report  describe  the  total  strain  increment  pro- 
cedure utilized  in  the  general  soil  cap  model.  A computer  code  was  written 
to  evaluate  the  material  constants  for  a prescribed  total  strain  path,  and 
a computer  subroutine  based  on  the  total  strain  increment  procedure  was  im- 
plemented in  a one-dimensional  wave-propagation  code.  Comparisons  of  re- 
sponses calculated  with  the  specified  cap  model  and  another  constitutive 
model  were  made  for  several  materials  subjected  to  shock  loadings.  Finally, 
the  limitations  and  problems  associated  with  the  cap  model  were  addressed 
and  recommendations  for  future  use  were  developed. 
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SECTION  II 
SOIL  CAP  MODEL 


TOTAL  STRAIN  INCREMENT  PROCEDURE 


To  analyze  any  plasticity  model,  it  is  beneficial  to  develop  a method  for  eval- 
uating the  material  behavior  for  a prescribed  stress  or  strain  path.  Since 
many  wave-propagation  codes  utilize  the  total  strain  increment  procedure,  this 
method  was  used  to  study  the  soil  cap  model.  The  following  derivations  are 
based  on  the  assumption  that  the  total  strain  increments  are  known.  However, 
before  the  total  strain  increment  approach  is  described,  some  background  in- 
formation must  be  presented. 

In  classical  plasticity  theory,  the  components  of  the  plastic  strain  rate  ten- 
sor can  be  expressed  as 

^ <-J 

where 

de^y  = plastic  strain  rate  or  increment 
d\  = scalar  function 
^ = plastic  potential  function 
o^y  = stress  components 

If  the  plastic  potential  function  is  assumed  to  be  identical  to  the  yield  func- 
tion, f,  an  associated  flow  rule  is  obtained;  i.e., 

deP;  = d,\  (2) 


The  following  expressions  describe  the  usual 
finitesimal  strains: 


e 

hi 


plasticity  relationships  for  in- 


= C 


ijkl 


^%l  - 


hi 


(3) 
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t 
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where 


^ijkl 


de 

de 


e 

kl 

kl 


= compliance  components 

= elastic  strain  rate  or  increment 
= total  strain  rate  or  increment 


The  consistency  equation  for  the  yield  function  can  be  expressed  as 


3f 

,eP. 


deP.  = 0 
<-1 


where  f is  assumed  to  be  a function  of  stress  and  plastic  strain, 
equations  (2),  (3),  and  (4)  results  in  the  following  expressions: 


(4) 


Combining 


(5) 


(6) 


In  summary,  if  the  total  strain  increments  are  known,  dX  can  be  computed  from 
equation  (6).  The  plastic  strain  increments  are  then  obtained  from  equation 
(2)  and  the  elastic  strain  increments  are  determined  by  subtracting  the  plas- 
tic strain  increments  from  the  total  strain  increments.  Finally,  the  stress 
increments  are  derived  from  equation  (3). 


YIELD  CRITERIA 

Figure  1 shows  the  general  soil  cap  model  yield  surface  in  two-dimensional 
stress  space.  The  elastic/plastic  model  is  rate-independent  and  combines  ideal 
plasticity  and  strain  hardening.  In  its  most  general  form,  the  yield  criterion 
is  a function  of  all  the  stress  components  and  the  stress/strain  history  of  the 
material . 
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For  the  particular  cap  model  studied,  the  ideally  plastic  portion  of  the  yield 
surface  was  assumed  to  be  of  the  form 

VjJ)  = 1 (7) 

where  J is  the  first  invariant  of  the  stress  tensor  and  j'  is  the  second  in- 
1 2 

variant  o^  the  stress  deviator  tensor.  For  the  cap  portion  of  the  yield  sur- 
face, the  expression 

= 1 (8) 

was  used.  Both  portions  of  the  yield  surface  were  used  in  conjunction  with  an 
associated  flow  rule  to  satisfy  Drucker's  stability  postulate  (ref.  4). 

Most  soil  behavior  can  be  adequately  described  by  using  yield  criteria  in  the  form 
of  equations  (7)  and  (8).  Weidlinger  Associates  specify  the  yield  criteria  in 
a slightly  different  form;  fj(Jj,  ""  0 ^nd  f2(Ji. 

represents  the  volumetric  plastic  strain  and  the  cap  expands  or  contracts  in 
stress  space  as  a function  of  eP.  Compaction  of  the  material  results  in  the 
expansion  of  the  cap  to  the  right.  Since  strain  hardening  can  be  reversed,  the 
soil  cap  model  allows  control  of  the  dilatancy  (inelastic  volume  change)  of  the 
material.  A generalized  form  of  Hooke's  Law  is  used  to  describe  the  elastic 
behavior  of  the  material  in  loading  and  unloading. 

The  ideally  plastic  portion  of  the  yield  surface  used  in  this  study  is  described 

]fK'>  = 

where  a,  y»  and  & are  material  constants.  This  exponential  form  of  the  yield 
surface  is  used  to  fit  the  data  obtained  from  triaxial  stress  experiments. 


4.  Drucker,  D.  C.,  "On  Uniqueness  in  the  Theory  of  Plasticity,"  Quav.  Appl. 
Math.  , 14,  1956,  pp.  35-42. 
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It  was  assumed  that  the  cap  portion  of  the  yield  surface  could  be  described  by 
a family  of  ellipses;  i.e.. 


f (J 
2 1 


1 


R'b' 


R^J 


(J, 


- C)' 


(10) 


where 


and 


Rb  = C - X (11) 

R = ratio  of  major  to  minor  axis  of  ellipse 
X = Jj  value  at  intersection  of  cap  and  J^-axis 
C = value  at  center  of  ellipse 
value  at  Jj  = C 

These  ellipses  have  horizontal  tangents  at  their  intersection  with  the  ideal 
yield  surface.  This  forces  the  plastic  strain  rate  vector,  which  is  always 
perpendicular  to  the  yield  surface,  to  be  vertical  at  the  intersection  and, 
thus, precludes  further  cap  motion  and  controls  the  amount  of  dilatancy. 

The  function  X depends  on  the  volumetric  plastic  strain  and  was  assumed  to  be 
of  the  form 


x(cP)  = 


£m 


D + Z 


(12) 


for  this  study.  Hydrostatic  test  data  for  a sand  were  used  by  DiMaggio  and 
Sandler  (ref.  5)  to  determine  the  specific  form  of  equation  (12),  but  other 
forms  may  be  appropriate  if  they  give  a better  representation  of  the  hydro- 
static behavior  of  the  material.  Additional  constants  that  can  be  evaluated 
from  various  laboratory  experiments  are  R,  D,  W,  and  Z.  Details  of  the  com- 
posite yield  surface  are  shown  in  figure  2. 


STRESS  AND  STRAIN  GRADIENTS  OF  YIELD  FUNCTIONS 

It  is  necessary  to  compute  the  stress  and  strain  gradients  of  the  yield  functions 
in  a total  strain  increment  procedure  before  equation  (6)  is  applied.  For  the 


5.  DiMaggio,  F.  L.,  and  Sandler,  I.  S.,  "Material  Model  for  Granular  Soils," 
Jourruxl  of  the  EnjincenCruj  Meahanias  Divieion,  Proceedings  of  the  American 
Society  of  Civil  Engineers,  Vol.  97,  No.  EM3,  June  1971,  pp.  935-949. 


ideally  plastic  portion  of  the  yield  surface,  only  the  stress  gradient  is  re- 
quired since  f^  is  not  a function  of  strain;  but,  both  stress  and  strain  gra- 
dients are  needed  for  the  cap  portion  of  the  yield  surface. 


The  stress  gradient  of  f^  can  be  expressed  as 


I 

J 3a. . 


VT’yeeW' 


- ^fy[.  - ve«i  j,  . 

Where  . are  the  components  of  the  deviatoric  stress  tensor  and  6^.y  is  the 
Kronecker  delta;  f^  is  defined  in  equation  (9).  Similarly,  the  stress 
gradient  of  fz  [equation  (10)]  is 

R'a*? . + 2(J  - C)6;  ; (1 

R^b^L  •'J 


Before  an  explicit  form  of  the  strain  gradient  of  f^  can  be  obtained,  some  in- 
termediate partial  derivatives  must  be  calculated.  Since  the  function  X [eq. 
(12)]  depends  solely  on  e^,  the  strain  gradient  of  X can  be  written  as 


apP  I eP 

DW  1 + ^ 


Equation  (9)  can  be  used  to  obtain  the  expression 


and  hence 


b = a - ye 


1 3C 


where  equation  (M)  has  been  used.  Rearranging  equations  (16)  and  (17)  gives 

_ax 

3C  _ 


9eP. 


aePy  (l  + Ryge®^) 


6.  . 
-LL 


DU 


eP 

1 + — 
' W 


1 + RySe 


6C 


(18) 


The  strain  gradient  of  b can  then  be  written  as 

3b 


Se?. 


ij  DW(l+§^ 


1 + RyBe 


BC 


(19) 


Combining  equations  (10),  (12),  (17),  and  (18)  gives  the  expression  for  the 
strain  gradient  of  fz  as 


- 2 


9eP. 


9Jl 


R'b^ 


pj;  * (0.  - o]  - ^ (J.  - c) 


- - ' ^^[f,R^bg(b  - g)  ^ (J,  - O] 

R'b'DW[l  + 6(Ra  + X - C)] 


9eP. 


(20) 


The  product  of  the  strain  and  stress  gradients  is  also  required  in  equation  (6) 
and  this  product  is 


9o . . 


where 


-scVX*  ■ 

[c  + 2f^6| 

1'  ■ f)] 

DW 

c = 

[l  + B(Ra  + X - o] 
2(Jj  - C) 

7 2 

R b 


(21) 


EVALUATION  OF  MATERIAL  CONSTANTS 


Experimental  data  from  triaxial  stress  experiments  are  required  to  evaluate 
the  material  constants  which  describe  the  yield  function  fj.  It  is  assumed 
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that  a judicious  choice  of  the  constants  a,  y.  and  3 will  define  a yield  sur- 
face which  will  adequately  match  the  experimental  yield  data. 

For  large  values  of  -J  , equation  (9)  becomes 


and  hence  a,  which  represents  the  von  Mises  Limit,  can  be  determined  if  the 
VoJ  values  of  the  experimental  data  approach  a constant  as  -J^  increases.  If 
Jj  = 0 the  yield  function  f^  is 

1 a - Y 

The  quantity  (a  - y)  is  proportional  to  the  cohesion  of  the  material,  and  y 
can  be  determined  once  a is  chosen.  The  constant  3 controls  the  degree  of 
curvature  of  the  yield  function  f^.  Hence,  3 can  be  chosen  to  give  the  best 
approximation  of  the  experimental  data  as  VjT  transitions  from  (a  - y)  to  a. 


A hydrostatic  test  is  required  to  evaluate  the  material  constants  which  de- 
scribe the  cap  portion  of  the  yield  surface.  If  a material  is  loaded  hydro- 
statically, the  stress  field  is  defined  by  normal  stress  components  that  are 
all  equal  to  the  negative  of  the  confining  pressure,  P,  and  shearing  stress 
components  that  are  identically  zero.  This  hydrostatic  behavior  is  repre- 
sented in  the  cap  model  as 


Equation  (12)  then  becomes 


X = -3P 


1 uP 

-3P  = Z + ^ fn  1 - ^ 


= W 


+ 3P)j 


where  yP  is  the  plastic  volumetric  compression.  The  total  volumetric  compres- 
sion, y,  is  then  expressed  as 


y = y®  + yP 


^*w[l 
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where  u®  is  the  elastic  volumetric  compression  and  K is  the  bulk  modulus.  The 
unloading  portion  of  the  hydrostat  is  assumed  to  be  elastic  and,  hence,  K is  de- 
termined from  the  slope  of  this  straight  line  (fig.  3).  Any  other  unloading 
test,  such  as  uniaxial  strain,  can  be  used  to  evaluate  the  other  required  elas- 
tic constants.  Evaluation  of  the  parameter  Z is  also  required  and  its  value 
can  be  estimated  if  a seismic  toe  exists  in  the  hydrostat  (fig.  3). 

The  parameters  D and  W determine  the  shape  of  the  loading  portion  of  the  hydro- 
stat. If  it  is  desired  to  match  the  point  on  the  hydrostat  (pi.  Pi),  D and  W 
must  satisfy  the  relationship 

Consequently,  D and  W can  be  varied  until  they  produce  a best  fit  of  the  load- 
ing portion  of  the  hydrostat.  Variations  in  D and  W change  the  degree  of  con- 
cavity of  the  loading  hydrostat.  Hence,  several  plots  of  p versus  P can  be 
constructed  for  the  various  choices  of  D and  W and  a final  engineering  judgment 
can  be  made. 

The  final  material  constant  (R)  is  determined  by  a trial -and-error  process. 

An  iterative  procedure  must  be  used  to  obtain  the  R value  which  gives  results 
that  are  compatible  with  experimental  data.  Values  of  R in  the  range  of  2.5 
to  5 have  been  used  and  have  produced  satisfactory  results. 


SECTION  III 


COMPUTATIONAL  PROCEDURES 


The  application  of  the  cap  model  to  wave-propagation  problems  was  a two-fold 
process.  Initially,  a computer  code  which  utilized  the  soil  cap  constitutive 
equations.  Program  ELLSTR  (appendix  A),  was  written  so  that  the  material  con- 
stants could  be  evaluated  for  a prescribed  total  strain  path.  An  equation-of- 
state  subroutine  based  on  the  cap  model  (appendix  B)  was  then  implemented  in  a 
one-dimensional,  wave-propagation  code.  Results  were  obtained  from  several  ma- 
terial models  subjected  to  simulated  shock  loadings  and  these  data  were  compared 
with  the  results  obtained  from  another  constitutive  model  and  actual  experimen- 
tal data. 

Program  ELLSTR  computes  stresses  and  elastic  and  plastic  strain  components  for 
a prescribed  total  strain  path.  This  program  is  based  on  the  procedure  outlined 
in  section  II  for  prescribed  total  strain  increments.  The  total  strain  increment 
procedure  for  the  soil  cap  constitutive  model  is  described  as  follows. 


The  material  constants  and  a parameter  that  controls  the  degree  of  subincremen- 
tation of  the  total  strains  are  required  as  input  to  ELLSTR.  The  user  also 
supplies  sets  of  total  strains,  and  the  difference  between  two  successive  sets 
of  total  strains  provides  the  strain  increment  set. 

Initially,  the  incremental  strains  are  assumed  to  be  elastic  and  the  stress 
components  are  computed  on  this  assumption.  Stresses  are  computed  from  the 
isotropic  stress/strain  relationship 

(a)  = [C]  {ee} 

where  {o}  and  {e®}  are  the  stress  and  elastic  strain  vectors,  respectively, 
and  [C]  is  the  constitutive  matrix.  Once  the  stresses  have  been  calculated 
the  appropriate  yield  function,  fj  or  f^,  is  evaluated.  If  the  calculated 
value  of  is  greater  than  C,  f^  is  evaluated  and  if  is  less  than  C,  f^, 
is  evaluated.  If  the  value  of  the  yield  function  is  less  than  1,  the  assump- 
tion of  elastic  behavior  is  presumed  to  be  correct  and  the  program  reads  in 
the  next  set  of  total  strains.  If  the  yield  function  is  greater  than  1,  some 
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plastic  deformations  must  occur  to  prevent  f from  exceeding  1.  For  this  case, 
the  total  strains  are  reset  to  their  previous  values  and  the  strain  increments 
are  then  subdivided  into  a number  of  smaller  increments.  The  number  of  sub- 
divisions is  proportional  to  the  value  of  the  computed  yield  function. 

For  any  particular  set  of  subincremental  strains,  the  stresses  and  yield  func- 
tion are  again  evaluated,  and  if  the  resulting  f is  less  than  1,  the  subincre- 
ment of  strains  is  presumed  to  be  elastic  and  the  program  proceeds  with  a new- 
set  of  total  strain  subincrements.  If  the  computed  f exceeds  1,  an  iteration 
loop  is  entered  in  which  new  stresses  and  plastic  strains  are  calculated  based 
on  equations  (1)  through  (6).  This  loop  terminates  when  f meets  some  criter- 
ion for  satisfying  the  yield  condition.  The  program  uses  a tolerance  limit  of 
0.99  < f 1.01  as  the  specified  criterion. 

There  is  a possibility  that  a value  of  f less  than  1 can  be  calculated  during 
the  iteration  loop  and  f can  oscillate  from  a value  greater  than  1 to  a value 
less  than  1.  Significant  errors  can  accumulate  if  the  iteration  loop  termi- 
nates with  a yield  function  value  less  than  1.  A method  has  been  adopted  in 
the  program  to  prevent  this  type  of  behavior.  It  is  noted  when  a value  of  f 
less  than  1 is  first  reached,  and  then  the  plastic  strain  increments  are  re- 
duced by  10  percent  successively  until  f reaches  1 from  outside  the  yield  sur- 
face. Actually,  if  f satisfies  the  tolerance  criterion,  the  iteration  loop  is 
terminated. 

After  each  calculation  of  the  plastic  strain  increments,  the  new  position  of 
the  cap  must  be  computed.  The  cap  parameter  X is  computed  from  equation  (12) 
but  the  transcendental  equation 

C = X + Ra  - 

must  be  solved  to  obtain  the  value  of  C.  A subroutine  which  incorporates 
Newton's  iterative  procedure  is  used  to  solve  this  transcendental  equation. 

The  other  cap  parameter  (b)  can  be  calculated  once  X and  C are  known. 

The  program  computes  new  elastic  and  plastic  strain  components  and  stresses 
for  each  set  of  total  strain  subincrements.  After  these  calculations  have 
been  completed  for  the  total  number  of  subincrements,  a new  set  of  total 
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Figure  4.  Flowchart  of  One  Cycle  of  Total  Strain 
Increment  Procedure  for  Soil  Cap  Model 
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strains  is  read  in  and  the  total  strain  increment  procedure  is  repeated.  A 
flowchart  describing  the  essential  features  of  the  total  strain  increment  pro- 
cedure for  a soil  cap  constitutive  model  is  shown  in  figure  4. 

The  soil  cap  model  was  implemented  in  the  one-dimensional,  wave-propagation 
code,  WONDY  IV  (ref.  6).  Since  this  code  utilizes  the  total  strain  increment 
approach,  the  essential  features  of  ELLSTR  were  directly  incorporated  in  the 
WONDY  IV  Code.  However,  modifications  were  made  in  the  coding  to  account  for 
the  absence  of  shear  strains  in  one-dimensional  motion. 


6.  Lawrence,  R.  J.,  and  Mason,  D.  S.,  WONDY  IV  ~ A Computer  Program  for  One- 
Dimensional  Wave  Propagation  with  Rezoning,  SC-RR-71 -0284,  Sandia  Labora- 
tories, Albuquerque,  New  Mexico,  August  1971. 
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SECTION  IV 
COMPARATIVE  RESULTS 


McCORMICK  RANCH  SAND 

Results  obtained  with  the  specific  cap  model  were  compared  with  the  behavior 

of  real  geologic  materials.  Since  data  were  available  for  McCormick  Ranch 

sand  (refs.  7 and  8)  Weidlinger  Associates  determined  the  cap  model  parameters 
for  this  material  from  the  uniaxial  strain,  triaxial  stress,  and  hydrostatic 
data.  DiMaggio  and  Sandler  (ref.  5)  used  these  data  to  estimate  the  following 
cap  parameters  for  McCormick  Ranch  sand: 

a = 250  psi  (1.72  MPa)  R = 2.5 

B = 6.7  X 10-Vpsi  (9.7  x lO'VPa)  D = 6.7  x 10-Vpsi  (9.7  x 10-VPa) 

Y = 180  psi  (1.24  MPa)  W = 0.066 

Young's  modulus,  E,  and  Poisson's  ratio,  v,  were  estimated  to  be  100  ksi  (689 
MPa)  and  0.25,  respectively.  The  laboratory  data  did  not  show  a noticeable 
elastic  region  upon  initial  loading;  therefore,  the  value  of  the  cap  parameter 
Z was  arbitrarily  chosen  to  be  -5  psi  (-0.034  MPa).  These  parameters  were  used 
in  ELLSTR  to  simulate  uniaxial  strain  and  hydrostatic  behavior. 

Figure  5 shows  a comparison  of  the  experimental  data  and  the  cap  model  results 
for  uniaxtal  strain  and  hydrostatic  tests.  These  cap  model  results  compare 
quite  well  with  the  laboratory  uniaxial  material  behavior,  although  there  is  a 
noticeable  difference  in  the  two  responses  in  the  unloading  portion  of  the  test 
The  laboratory  data  exhibit  a more  predominant  recovery  of  strain  than  do  the 
cap  model  results.  Weidlinger  Associates  reported  that  the  hydrostatic  data 
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7.  Zelasko,  J.  S. , and  Ingram,  J.  K. , Soi.t  Property  Investigation  and  Free- 
Field  Ground  Motion  Measurements,  Project  BACKFILL,  USAEWES,  Vicksburg, 
Mississippi,  December  1967. 

8.  Mazanti , B.  B. , and  Holland,  C.  N.,  Study  of  Soil  Behavior  Under  High  Pres- 
sure, Report  1:  Response  of  Two  Recompacted  Soils  to  Various  States  of 
Stress,  Report  S-70-2,  USAEWES,  Vicksburg,  Mississippi,  February  1970. 
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for  McCormick  Ranch  sand  were  unreliable  because  of  possible  errors  in  the 
radial  strain  measurements  (ref.  5).  Nevertheless,  the  cap  model  described 
in  this  report  was  used  to  simulate  hydrostatic  behavior  and  these  results 
are  also  shown  in  figure  5. 

A CIST  experiment  was  simulated  to  evaluate  the  response  of  the  cap  model  in  a 
field  loading  condition.  Another  constitutive  model,  one  in  which  a piecewise 
linear,  hysteretic  hydrostat  and  an  isotropic  yield  surface  are  used  to  describe 
the  material  behavior,  was  also  used  in  the  CIST  simulation.  This  model,  here- 
after referred  to  as  the  engineeriruj  models  has  been  used  in  the  analyses  of 
several  CIST  events.  Reference  9 contains  a complete  description  of  this  model 
and  a typical  CIST  event. 

Piecewise  linear  representations  of  the  hydrostat  and  yield  surface  for 
McCormick  Ranch  sand  were  needed  for  the  engineering  model.  Therefore,  the 
hydrostat  and  yield  surface  obtained  from  the  cap  model  were  approximated  by 
straight-line  segments.  The  cap  model  and  engineering  model  representations 
of  the  McCormick  Ranch  sand  hydrostat  and  yield  surface  are  shown  in  figure  6. 

The  WONDY  IV  Code  was  used  in  the  simulation  of  a CIST  event  at  a relatively 
shallow  depth.  A typical  CIST  cavity  with  a 12-inch  (30.5-cm)  radius  and  an 
assumed  pressure  function  of 

P = 5800  psi  (40  MPa)e‘^^°^ 

where  P = pressure  and  t = time  (sec)  were  used  in  the  analysis.  The  materials 
were  modeled  in  cylindrical  geometry  with  variable  zone  sizes.  Zones  closest 
to  the  explosive  cavity  were  approximately  3 inches  (7.6  cm)  wide  and  the  zones  ex- 
panded to  a maximum  width  of  18  inches  (45.7  cm)  at  a range  of  125  feet  (38.1  m). 
This  mathematical  simulation  of  the  CIST  event  was  performed  for  a duration  of 
20  msec. 

Time  histories  of  the  particle  velocities,  radial  stresses,  and  other  responses 
were  obtained  in  the  CIST  simulation  with  both  the  cap  and  engineering  models. 

9.  Bratton,  J.  L.,  Fedock,  J.,  and  Higgins,  C.  J.,  A Pavametria  Study  of  the 
Effects  of  Material  Properties  Upon  Cylindrical  Wave  Proixigation , AFWL 
Technical  Note  (in  preparation). 


22 


Figure  7 shows  the  velocity/time  histories  obtained  at  four  ranges.  Compari- 
sons of  the  radial  stress/time  histories  calculated  with  both  models  are  shown 
in  figure  8,  There  is  very  little  difference  in  the  responses  obtained  with 
the  cap  and  engineering  models.  The  arrival  times  and  loading  phases  of  the 
time  histories  computed  with  both  models  show  no  appreciable  differences.  At 
all  four  ranges  the  peak  velocities  calculated  with  the  cap  model  are  slightly 
higher  than  the  amplitudes  obtained  with  the  engineering  model.  However,  the 
velocity/time  histories  computed  with  the  cap  model  decay  more  rapidly  than 
the  corresponding  responses  from  the  engineering  model.  Radial  stresses  cal- 
culated with  both  models  compare  quite  well,  although  at  the  more  distant 
ranges  the  cap  model  results  show  slightly  higher  stresses  than  do  the  engi- 
neering model  results.  Comparison  of  tangential  stresses  obtained  with  both 
models  shows  the  same  phenomenon  as  did  the  radial  stress  comparison. 

STIFF  CLAY 

To  evaluate  the  adequacy  of  the  cap  model  in  estimating  in-situ  material  prop- 
erties, data  from  a typical  CIST  event  were  analyzed  with  the  cap  model.  An 
iterative  procedure  was  used  to  estimate  the  cap  model  parameter  values  that 
best  describe  the  dynamic  material  behavior.  These  parameters  were  then 
used  to  calculate  the  material  response,  and  the  results  were  compared  with 
experimental  records.  The  iterative  process  was  repeated  until  the  calculated 
responses  were  in  good  agreement  with  the  experimental  data.  Comparisons  of 
the  recorded  data  and  the  responses  computed  with  the  i«-s£tu-based  engineering 
model  were  also  made. 

One  geologic  layer  (stiff  clay)  at  a CIST  site  was  analyzed  with  the  cap  model. 
The  in-situ  material  model  parameters  were  derived  from  the  horizontal  velocity 
records  for  a depth  of  5 feet  (1.52  m). 

in-situ  values  of  the  density  and  Poisson's  ratio,  which  were  used  in  the  cap 
model,  were  based  on  laboratory  estimates  of  these  properties.  The  bulk  modu- 
lus was  estimated  from  the  unloading  portion  of  the  hydrostat  obtained  from 
laboratory  data.  Initial  estimates  of  cap  parameters  D and  VI  were  obtained 
by  determining  a best  fit  to  the  loading  portion  of  the  laboratory  hydrostat. 
The  value  of  R was  chosen  to  be  4.0.  Since  laboratory  estimates  were  available 
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Figure  8.  Comparison  of  Radial  Stress/Time  Histories  from  Cap  and  Engineering  ? 

Model  Representations  of  McCormick  Ranch  Sand  f 
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for  the  yield  surface,  values  of  a,  B,  and  y were  initially  chosen  to  fit 
these  yield  data. 


Once  initial  estimates  of  the  cap  model  parameters  were  established,  the  WONDY 
IV  Code  was  used  to  compute  the  material  response  based  on  an  assumed  input 
pressure  function.  After  the  cap  model  parameters  were  varied  several  times 
during  the  iterative  process,  the  parameters  which  produced  the  best  fit  to 
the  experimental  velocities  were  obtained.  These  were  as  follows: 

a = 46  psi  (0.32  MPa)  0 = 4.0  x lO'Vpsi  (5.8  x 10-VPa) 

B = 2.5  X 10-Vpsi  (3.6  x lO'VPa)  W = 0.040 

Y = 12  psi  (0.08  MPa)  v = 0.32 

R = 4.0  K = 2.6  X 10®  psi  (1790  MPa) 

The  responses  obtained  with  the  initial  estimates  of  the  cap  model  parameters  [ 

compared  fairly  well  with  the  experimental  data;  i.e.,  the  final  estimates  of 
many  parameters  were  not  significantly  different  from  their  initial  values. 

In  particular,  only  D and  W,  which  were  initially  estimated  to  be  1.0  x 10”V 
psi  (1.45  X 10"’/Pa)  and  0.031,  respectively,  were  greatly  varied  in  the  iter-  f 

ative  process.  A similar  iterative  process  was  used  to  estimate  the  parameters  ! 

for  the  engineering  model,  which  was  also  used  to  describe  the  in-situ  behavior  | 

of  the  stiff  clay.  I! 

Tensile  properties  of  the  stiff  clay  could  not  be  absolutely  determined  be-  i 

cause,  in  general,  the  CIST  event  did  not  produce  material  responses  with  li 

noticeable  tensile  phases.  Hence,  no  tensile  properties  of  the  stiff  clay  | 

were  estimated.  l 

In  both  iterative  processes  the  assumed  input  pressure  function  was  ! 

P = 5800  psi  (40  MPa)e'^®°^ 

The  geometric  modeling  of  the  material  for  both  cases  was  the  same  as  that 
described  previously  for  the  McCormick  Ranch  sand.  A simulation  time  of  30 
msec  was  used  for  both. 


27 

'"''-''i”'..  A 


Final  estimates  of  the  cap  model  parameters,  which  described  the  clay  mate 


rial  at  the  5-foot  (1.52-m)  depth,  were  used  as  input  to  ELLSTR.  Hydrostatic 
behavior  was  simulated  in  the  program  and  the  results  were  compared  with  the 
hydrostats  obtained  with  the  in-situ  engineering  model  and  the  laboratory- 
inferred  properties  (fig.  9),  The  cap  model  parameters  (a,  3,  and  y)  were 
used  to  describe  the  final  estimate  of  the  yield  surface.  A comparison  of 
the  estimated  in-situ  yield  surfaces  obtained  with  the  cap  and  engineering 
models  and  the  laboratory-based  properties  is  also  shown  in  figure  9. 

Velocity/time  histories  calculated  with  the  cap  and  in-situ-based  engineering 
models  are  compared  with  the  experimental  records  in  figures  10  and  11,  re- 
spectively. A comparison  between  the  responses  calculated  with  the  laboratory- 
based  engineering  model  and  the  field  records  is  shown  in  figure  12. 

Arrival  times  calculated  with  both  in-situ  models  are  in  good  agreement  with 
the  field  data.  Peak  velocities  computed  with  the  cap  model  are  lower  than 
the  corresponding  amplitudes  calculated  with  the  engineering  model.  The 
loading  phases  and  peak  velocities  calculated  with  the  cap  model  match  the 
experimental  data  quite  well  at  the  3-foot  (0.91-m)  and  5-foot  (1.52-m)  ranges, 
but  the  computed  amplitude  is  approximately  40  percent  higher  than  the  re- 
corded amplitude  at  the  8-foot  (2.42-ni)  range.  Neither  in-situ  model  gives  a 
good  representation  of  the  actual  material  behavior  at  the  8-foot  (2.42-m) 
range.  The  unloading  portions  of  the  time  histories  calculated  with  both 
in-situ  models  compare  reasonably  well  with  the  experimental  data.  In  gen- 
eral, the  differences  between  the  responses  calculated  with  both  in-situ 
models  are  not  significant. 
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Figure  9.  Comparison  of  Cap  and  Engineering  Model  Representations 
of  In-Situ  Properties  and  Laboratory- Inferred  Properties 
for  Stiff  Clay 


Pressure,  MPa 


Velocity,  ft/sec  Velocity,  ft/sec  Velocity,  ft/sec  Velocity,  ft/sec 


P 


j 

Figure  10.  Comparison  of  Velocity/Time  Histories  from  Estimated  in-Situ 

Cap  Model  Parameters  and  CIST  Data  for  Stiff  Clay  1 

30  ] 


4 


Velocity,  m/sec  Velocity,  m/sec  Velocity,  m/sec  Velocity,  m/sec 


Velocity,  ft/sec  Velocity,  ft/sec  Velocity,  ft/sec  Velocity,  ft/sec 


P 


Figure  11.  Comparison  of  Velocity/Time  Histories  from  Estimated  in-Situ 
Engineering  Model  Parameters  and  CIST  Data  for  Stiff  Clay 
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Figure  12.  Comparison  of  Velocity/Time  Histories  from  Laboratory-Based 
Engineering  Model  Parameters  and  CIST  Data  for  Stiff  Clay 
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SECTION  V 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  cap  model  is  capable  of  representing  laboratory  and  dynamic  in-situ  behavior 
for  at  least  two  distinct  geologic  materials.  Results  from  an  analysis  of  a 
CIST  event  show  that  the  cap  model  can  be  used  to  produce  velocity/time  his- 
tories that  are  in  reasonably  good  agreement  with  experimental  data. 

The  cap  model  meets  all  uniqueness  and  stability  requirements.  (The  theoreti- 
cal basis  for  the  cap  model  is  described  in  reference  5.)  Since  the  cap  model 
utilizes  an  associated  flow  rule  with  a convex  yield  surface,  it  satisfies 
Drucker's  stability  postulate;  the  engineering  model  uses  a nonassociated  flow 
rule  and  does  not  satisfy  the  postulate. 

The  cap  model  described  in  this  report  and  implemented  in  the  WONDY  IV  Code  has 
several  disadvantages,  however.  The  procedure  for  fitting  both  laboratory  and 
field  data  is  fairly  complex.  A trial -and-error  process  is  required  to  deter- 
mine values  of  the  cap  model  parameters,  but  recent  studies  (refs.  10  and  11)  have 
indicated  that  the  iterative  procedure  can  be  automated.  It  is  also  difficult 
to  relate  the  cap  model  parameters  to  the  material  properties  in  commonly  under- 
stood terms.  The  cap  model  equation-of-state  is  more  complex  than  the  equation- 
of-state  based  on  the  engineering  model.  However,  the  required  execution  times 
for  the  calculations  performed  in  this  study  were  not  appreciably  different  for 
the  two  models. 

Some  aspects  of  typical  CIST  responses  are  difficult  to  reproduce  with  the  cap  mod- 
el, but  this  difficulty  is  common  to  the  engineering  model.  The  CIST  data  often 


10.  Isenberg,  J.,  Collins,  J.  D.,  and  Kennedy,  B.,  Statistical  Estimations  of 
Geological  Material  Model  Parameters  from  Cylindrical  In-Situ  Test  Data, 

AFWL-TR-76-187,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force  Base, 

New  Mexico. 

11.  Isenberg,  J.,  Automated  Data  Analysis  Application  of  Statistical  Estima- 
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show  velocities  that  maintain  high  amplitudes  for  relatively  long  durations, 
and  both  the  cap  and  engineering  models  have  difficulty  in  matching  this  type 
of  behavior.  Additionally,  many  CIST  velocity  records  show  a sharp  decrease 
in  amplitude  after  the  initial  peak,  followed  by  a gradual  increase  in  the  ve- 
locity over  several  milliseconds.  This  behavior  is  also  difficult  to  match 
with  either  the  cap  or  engineering  model.  The  tensile  behavior  of  soils  is  an 
area  of  much  uncertainty,  and  neither  the  cap  nor  the  engineering  model  appears 
to  give  a completely  adequate  representation  of  this  type  of  behavior. 

For  the  cases  considered  in  this  study,  it  appears  that  if  the  cap  model  param- 
eters were  chosen  to  match  the  hydrostat  and  yield  surface  of  the  engineering 
model,  the  differences  between  the  time  histories  calculated  with  both  models 
would  be  minimal.  This  suggests  that  other  approaches  to  the  modeling  of  dy- 
namic soil  behavior  need  to  be  considered.  The  use  of  a residual  yield  surface 
may  be  one  means  for  obtaining  the  desired  waveform  behavior. 

Although  the  cap  model  satisfies  all  theoretical  requirements  and  can  adequate- 
ly represent  most  types  of  static  and  dynamic  soil  behavior,  it  does  not  appear 
that  the  cap  model  offers  significant  advantages  over  the  engineering  model  for 
the  cylindrical  cases  evaluated.  However,  the  guarantee  of  solution  uniqueness 
that  is  ensured  by  the  cap  model  may  provide  the  user  with  increased  confidence 
when  problems  involving  general  three-dimensional  conditions  are  calculated. 


APPENDIX  A 
PROGRAM  ELLSTR 


INPUT  REQUIREMENTS 

Card  1 - Format  (8A10) 

Title  Card  (up  to  80  characters) 

Card  2 - Format  (10F8.2) 

Material  Constants  - E,  v,  0,  W,  R,  a,  6,  y,  1 
Subincrement  Control  - PINC 

Card  3 to  Last  Card  - Format  (110,  6F10.2)  - NSTEP,ET(I) 

NSTEP  is  an  integer  to  identify  each  set  of  total  strains  that  is  read 
in.  Normally,  the  first  card  of  this  set  will  have  NSTEP  = 0,  and 
ET(I)  =0.,  0.,  0.,  0.,  0.,  0.  Subsequent  values  of  total  strain  define 
the  prescribed  total  strain  path. 

NSTEP  = -1  indicates  that  this  set  is  complete  and  that 
another  problem  set  is  next  in  line. 

NSTEP  = -2  indicates  that  this  set  is  the  last  complete 
problem  set. 

ET(I)  = components  of  total  strain. 


PROr.RAM  KLbTRt  iNPUt  tOUTPUl  vTAPE^*  iNPot  • I APEb«OOlPuT  ) 

COMMON  /blRkbb/blO(bl  tbDl  • bO?  • bU:?  tF  J I • F J2P  • F J2PH 
common  /SiRAIN/tPPL (6  » • tt L ‘ b • * t T < 6 * tOt T ( 6 > tOt PL < 6 * # tP I 
COMMON  /MATtR/C2  •C3«()2«CC  (6»6>  •ALPH*HLT  »oAM«K»Ui* 
common  /CAP/X*C  •UW.KiO*RB*Kbbg*RALPH*ROAM»C.Mttt1 
COMMON  /CONIRL^NCOUNT • XZ  * CH I bO iF Y •NO^M 
r 

OIKUNSION  T ! TlM  1C  » »F  TnE  W(6I 

r 

c this  PROORAm  AbbUMtb  A blRAlN  iNCRiMLNl  PRbSCRlBtO  PAlH 

C ALl  SIX  iNOtPLNOtNT  COMPONtNTS  QF  TOTAL  blRAlN  ARL  KtUUlRtO 


TO  RFAr>(  A*1000)  * M TlF  ( 1 » t I*1  tft) 

1000  FORmaT(8A101 

RFAn(  S.lOlOlf  *FnU*0*W*R*AI  PHtbF  T .CiAMtXlNl  *PtNC 
1010  FORMAT! 10F8*2) 

»RI  TE  <ft*  nOOM  I I iLf  U » *1  •!  *81 
uno  FORMAT!  lHl.l0X*aA10/  1 

«RlTF!6*inOU  «FNtl*D*M  •R*ALPH*HE  T iGAMvX  InI  tpl  nC 
1110  FORMAT! IH  •3Xi2Ht«»L8*2*3X*iHNO«*FA«2*3X»ZHU*»tB«2»iA*ZHN-»F6«A* 

1 3X«2HR»  tFb.Z*  3X«6HAlPHA* (Fb* I • 3X • bHbt T A « • L 8 • 2 • 3 A •bHGAMMA* *F b • i • 

2 3X,bHX!NT«  •F*>.  1 tlX.bHPiNC*  *rb,  1/  1 
C 

C INITIALUL  variaulls 

c 

CHISQ* (O.R/HFT ) •AlOG! ALPH/GAM ) 

FK«  ALPH-r,AM 

nw«f>»\N 

RSO«R»R 

ralph*r»alph 

RGAM«R*GAM 

gmbft«gam*bet 
X-XtnT 
XRAaX^RALPH 
C«0. 1 

CALL  COMPC!C*XRA*Rr,AM*RET  I 

RBT-X 

RR5iO*RR*RR 

XSQ*X*X 

G?«  F/!1.^FNU) 

C1-G2/I l.-2.*FNUl 
C2-C1*!  U-FNU) 

C3«r 1 •FNU 
FPI*0. 

FV»FJ1«FJ2PR»0, 

N INr.o 

ncount«o 
no  n t«i *6 
siGi  n«o. 

rT!  ! ) «o« 

F^L ! T )*0. 

FPPL ! ! »*0. 

no  lA  j*i*b 

IS  CC!  I 

no  ?0  !•! *3 
no  ?o  J*i *3 

20  CC! I • J)«C3 
no  72  1-1*3 
77  rr! T * I 1 -C2 
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no  24  T«4*A 

24  CC!  I • I )-C,2 

c 

C strain  iNCRtMtNl  CYCLt*  READ  NER  TqTAL  STRAINS 

r 

so  REAn ( S* lOSO INSTEP* It TNER! 11*1-1*61 
lOSO  FORMAT! 1 10*6F10.2» 

!F!NSTEP*E0«01G0  TO  700 
IFlNSTEP.EQ.-l  )CiO  TO  10 
|F!nSTFP,F0,-2)ST0P 
NCOONT-0 
r 

C COMPUTE  total  strain  INCREMENTS*  ASSUME  ELASTIC  ANU  CQMPulE  STRESSES* 

C 

no  SO  1-1*6 
OETl I )*ETNEW! I »-ET!  ! 1 
SO  FT! T I.ETNFW! I ) 

C CMECX  FOR  YIELDING*  IF!FY*LE*ll  INCREMENT  IS  ELASTIC 

C 


JL 
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NrNr*0 

NnUM«0 


CALL  YIFLO 
NnUM»l 

tF<FY.LE.U)GO  TO  700 
r 

C PFOUCF  SIZE  OF  TOTAL  STRAIN  INCREMENTS 

C 

ninc«ify-i. |#PINC 
I F I N I NC • G I • ZOOO ) N I NL  *2000 
|F<NINC*LT*1 )NINC-l 
FlNr*NlNC 
DO  100  I*l*6 
ETM  »aET<  n-OET(  I I 
100  DETI I IsDETI I )/FINC 

r 

C PROCEED  with  ThE  SHALL  TOTAL  STRAIN  INCREMENTS 
C 

no  ISO  K*l*NINC 
no  no  i«i*6 
no  ETC  I ).ET(  1 )♦DETC  I ) 

CALL  YIELD 

IFIFY.LE. 1.0)00  To  150 
CALL  PSTRN 
NINNER«0 
m CALL  YIELD 
r 

C 0.99  .LE.  FY  .lE.  l.Ol  IMPLIES  SATISFACTORY  STRESS  POINT 

IFlFY.GT.l.Ol )G0  TO  127 
1F<FY,GT.0.99)  GO  TO  ISO 
N!NNEf?*NlNNER^l 
IFCNINNER.GT.100)  GO  TO  ISO 
DO  l?ft  1*1.6 
FACT*-0.1 

EPPLC I )»EPPL* I )>FACT»OEPLC  I ) 

126  DEPLC l)«.9*0EPL( ll 

ePl»EPPLC 1 )^EPPL<2)^EPPLC3) 

X.XINT^C ALOGC l.O^EPl/W) )/n 
XRAsX^RALPH 


call  compccc.xra.rgam.pft ) 

RB*C-X 
RBSQ»RB*RB 
GO  TO  113 

127  NINNER*NINNER4-1 
lAO  CALL  PSTRN 
GO  TO  113 
C 

ISO  CONTINUE 
C 

730  WRITE! 6.1 150)NSTEP.NINC.NCOUNT 
mo  FORMATClM  • / . SX  .6HNSTEP*  » I 5 • SX  . 5HN  I NC*  • I 5 • SX  • 7MNC0UNT  ■ .15) 

r 

FMU*ETCI)4ET<2»4-ET(3) 

WRITE(6*1  UO)  CETC  I ) • I«i  .6)  •FMU 
1140  FORmATCIH  .4X.12HT0TAL  S I RA iN vAX . 6E 10. 3 . SX * 3HHU* .E 1 0. 3 ) 
r 

WRl TEC6. 1130) CEELC 1 ) . I «1 .6) 

1130  FORMATCIH  .4X.14HELAST IC  S TRA I N .2X .6E1 0. 3 ) 

r 

WRITEC6.  U2S)  CEPPL  C I ) . 1-1  .6)  .EPl 
112S  FORMATCIH  .4X * 1 4HPl AST IC  ST RA I N *2X «6E 10. 3 t SX »4HEP1 « .E iU. 3 ) 

C 

WRITEC6.1 120) CSIOC  J ) . 1 -1 .6) »FJ1 •FJ2PR 
1120  FORMAT! IH  . 1 OX . 6HSTRESS *4X tSE 10. 3 . SX . 3HJ 1 « * E U. S » SX »6HF J2PR«  tElA.S 
1 1 
C 

R-SORTCRBSO)/R 

WRITE! 6*1 122 )FK.B.C*X.FY 

1122  FORMAT! IH  . 4X . 2HY- .E 1 0. 3 • SA *2HB* . E 10. 3 . SX .2HL« .E 10. 3 tSA* 2HX- . 

1 E10.3*5X.3HFY«*E10.3) 

C 

GO  TO  SO 

STOP 

FNH 
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c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 


c 

c 

c 


c 


c 

c 

r 

r 

r 

r 


r 

r 

r 


r 


r 


SUAftOUTiNF  PSTRN 

COMMON  /STRESS/blG(6) 9^01 *b02*S03«F Jl «Fj2P«f J2PR 
COMMON  /STRA|n/£PPL(6»  tFEt  (6)  *ET(6>  •0en6)*0tPL(6»  ttPl 
common  /MATER/C2tC3»C>2*CC(6*6)  tALPHtUE T «GAM«R«0« « 
COMMON  /CAP/X»C*DW*RSO*RB*RBSU«RALPH*ROAM«GMaEt 

common  /contrl/ncount  tx^fCHiso.Fy indom 

OIMFNSION  GF(6) 


GF»GRA0IENT  of  F 
GF5»GRADIENT  of 

NCOUNT-NCOUNT+1 

IF«FJ1.GT,C»  GO 

yielding  on  cap 

COMPUTE  GF 


WRr  STRESS 
F MRT  strain 

To  20 

(F2I 


Al»FJl-C 

A?2*2,«Al 

A23«RSO/RBSO 

A24*A22/RRS0 


A25»2,»A23 

GFf 1 l■A21*Snl♦A24 

r>Fl2)»A23*Sn24A24 

GFI3)*A23»SD3>A24 

GF(4)«A25*SIGU) 

r,F(*i)«A?S»SrG(S) 

GF<AI«A2S»SIG(6» 

compute  GFS 

A12-X.C 

A3«A14FY»RB*8ET»IRB-RALPHI 
A4*RnS0«nW»l  1 •«-BET«(RALPH-fA12I  ) 
GFSS*-?«»EXP(-n*IX-XlNT» )«A3/A4 

R3«n* 

00  10  I « 1 • 3 
10  R3»R34^GFU)*GFSS 
r»0  TO  30 

YIELDING  ON  FI 

?0 


COMPUTE  GF 
AS«FY/I2**FJ2P) 

A6*2**AS»FY»GMBET«FJ2PR*CXPIBET«FJ1 } 

GF( 1 »»AS»SDl4A6 

GF(21«AS«S02^A6 

GFI1I»A5*SD3+A6 

A7«2««AS 

0F(4I«A7»S|G(4) 

GF<A1«A7*SIG(S1 

GF(6)>A7*SIGf6) 

COMPUTE  OLAM 

Pl»o. 

R2»n* 

DO  40  I«l»3 
00  40  J>1«3 

R1*B14CCI I *J)*0ET( I 1*GF( J) 

40  R?«B?4CCn  •J)*GFI  I I«GF(  Jl 

no  AO  I«4*A 

Bl*8142*«CC(  I • n»OET(  t l♦GF(l  \ 

AO  R2«B2^CC  I I # M»GFU  »«GF(  n 

0LAM«B1/(B2-B3I 


1 


1 

1 


I 

I 


i 

1 
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f 


( 

I 

t 

I 


TT 


C COMPUTE  PLASTIC  STRAIN 

r 

DO  60  t>l *6 

DEPL ( 1 I*DLAM*GFI I ) 

60  EPPlU  )»EPPLt  n^OEPL  ( I ) 

EPl»EPPL*  I )^EPPL*2I^EPPLI  -il 

r 

XsXINT4-< ALOGI  1«0>CP1/W)  I/D 
XRA>X4>RALPH 


CALL  COMPC(C*XRAfRGAM»BET) 

RB*C-X 

RRSO*Pfl*PB 

return 

END 


SUBROUTINE  YIELD 

common  /STRESS/SIGI6I tSOI •S02*S03*Fjl »FJ2P»FJ2PR 
common  /SIRAIN/EPPL (6) •ELL(6I vEI (6) »DET  <6) •DbPLl6) •LPl 
common  /MA1ER/C2tC3*02*CC(6*6l tALPHvBE T « OAm »R •U*n 
COMMON  /CAP/x#C#OW.RSO»RB*RbSO*RALPH»RGAM*&MbET 
common  /CONTRL/NCOUNT  tXZ  *CH1S0*FY  *NDUM 
C 

C COMPUTE  ELASTIC  STRAINS  AND  STRESSES* 

C 


DO  10  1*1 >6 

10  fCL( I l»ETI I |-EPPL< I I 

Sir,(  1 )»C2*EEL(1  )♦C3•^EFL<?I♦EEL(3)  ) 
SIGC  ?l  «C2»EEL<2I^C3«<EEL<  lUEEH  1 1 ) 
SlGm«C2*etL(3»^Cl*(EELM  ^♦EtL<2^ 
SIG(A»«G2«EELU» 

SIG(5)»G2*EEL(SI 

SIG(6l*G2»EeL(6l 


EJl »StG(l l+SIGI? I♦SIG<  31 
rj3*FJl/3. 

SD1*SIG( 1 I-FJ3 


Sn?«SlG(?>-FJ3 

SD3«StGl3)-FJ3 

FJ2P*CS01*S01>SD2»SD2+S03«SD3> /2.^ 

1 (S1G(4)**2)^<S1G(SI**2>'^(SIG<6I««21 
F J2PR*S0RT IFJ2P) 


C 

C 


C 

C 

20 


'»0 

1000 


1F(FJUGT*CI  GO  TO  20 


AaFJl-C 

FY*(RS0»FJ?P4.A*A)/RBS0 


IFInDUM*EO*OI  fy*sortifyi 


RETURN 

IFIFJI .GT.CHISOIGO  Tq  30 

DEN»ALPH-GAM*EXPIBET*FJ1 I 

FY*<Fj2PR1/DEN 

RETURN 

print  1000 

FORMAT! 10X*2AMWARnING  ON  FI  FROM  YIELD 
FJ1«0. 


FY«100. 

RETURN 


END 

subroutine  C0MPC!C*C1 *C?*C3) 
C 

C USE  NEWTON'S  METHOD 

C 

COLD-C 

NC«rt 

70  A1«C2»EXP!C3»C I 

A7«ri.c-Al 


I 


A3ar3»Al^U 

r«r«A2/A3 

IFIARS<C3*IC0L0-C) )*lT*0*0001 >00  TO  90 
NC«NC^l 

IE(NC*CiT«20IGO  TO  60 

COLn«C 

GO  TO  70 

10  WRITE(6*1100IC 

1130  EORMATIIH  tAX^lSNC  NOT  CONVERGING  iEJ0*3l 
STOP 

90  ie(c*gt.o*oic«o. 
return 

END 
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APPENDIX  B 

CAP  MODEL  EQUATION-OF-STATE  SUBROUTINE 

INPUT  REQUIREMENTS 

The  input  cards  are  essentially  identical  to  those  outlined  in  the  WONDY  manual. 
On  Card  2 set  NVAR  = 28,  and  on  Card  10  specify  2 for  the  equation-of-state 
Cards  15  and  16  contain  the  equation-of-state  constants  and  they  art  as  follow  ; 
Card  15  - Format  (7E10) 

CES  (1,  PLATE) CES(7, PLATE)  = initial  mass  density, 

initial  bulk  sound  speed.  Young's  modulus,  Poisson's 
ratio,  alpha,  beta,  and  gamma,  respectively. 

Card  16  - Format  (7E10) 

CES  (8,  PLATE) CES(10, PLATE)  = R,  D,  and  W,  respectively. 

Note:  The  coordinate  of  the  initial  cap  has  been  set  to  .02a 
in  the  program  and  subincrement  control  has  been  set  to  20. 
Appropriate  changes  can  be  made  in  the  coding  to  modify 
these  set  values. 

In  Subroutine  GENERAT  the  values  of  C,  X,  and  b must  be  preset  by  a call  to 
STIN2  (entry  point  in  STAT2)  and  stored  in  the  appropriate  locations  in  the 
vector  STORE. 

Subroutine  STAT2  predicts  stresses  and  strains  based  on  the  cap  model  equation- 
of-state.  In  the  first  segment  of  the  subroutine,  the  basic  material  constants 
are  set  up  if  a new  material  zone  is  involved,  and  then  the  previous  values  of 
stresses,  strains,  and  cap  coordinate  data  are  recalled. 

The  standard  WONDY  procedure  is  followed  in  computing  the  strain  rates,  with 
the  total  strain  increments  and  total  strains  evaluated  next.  Then  the  sub- 
routine closely  follows  ELLSTR  with  the  exception  that  the  evaluations  per- 
formed by  YIELD  and  PSTRN  are  incorporated  directly  into  STAT2. 

New  values  of  stresses,  strains,  cap  coordinates,  and  certain  other  variables 
are  saved  in  the  vector  DATB  before  the  return  from  the  subroutine  is  made. 
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subroutine  STAT2 
C WONnV  IVt  BARON  CAP  MOnfL 

COMMON  /CONST/  AODATAI 1001 ibl •B2>CESI8A>20I itXI T> lALPHAt 

1 lOUMPt  J«  JA>  KM(3|<  KTI>  L>  UMAX i COL i LUK>  LPHA> 

2 NOMESHt20)»  NONE*  NOP*  NTwO*  NVAR • PLATE*  PRINiE*  PRInTL* 

3 PRINTN*  STATEI20I*  SUMQE * TOEP*  WL * WR*  W402U*  X2EK0 

COMMON  /names/  C*  CN*  t*  £N*  M.  P.  PA*  PN*  0*  QA.  ON*  R*  RA*  RN* 

1 S*  SA*  SN*  U*  UA*  Un*  UBN*  UN*  X*  XA*  XB * XBN*  XN*  2i  LA*  2N 
COMMON  /retain/  N*  T.  OElTISI*  ETOT*  HToI*  IlAbLtlSOt*  IAbLEl2*50l 
1 * PERACTI I80SI * OERACT ( laOSI * LACT 

COMMON  /JMSO/  DATbUOOI*  OELRMO*  DELXJ*  JONE  * NEkPLAT.  RHODOT 
COMMON  /»JSO/  Dtp.  ECONST(20I*  FCONSTIIpOI.  ECK1II2UI*  ECRITitZO). 

1 JIAPL*  LDUMP*  NSTARI*  SIOACT*  SIGMAF(;0I*  SIOMAIE I2UI • 

2 SIGMA0(20)*  SIGMAOI(20)*  SIGSEP 
TYPE  INTEGER  PLATE*  NA020 

TYPE  REAL  KM.  kTI.  M,  NOMESH 

r 

C SET  UP  variables  FOR  NEW  PLATE 

f 

IFiNFMPLATlS.IO 
S RHOn«CES( 1 .plate  I 
fO*rFS(2*PLATEI 
H'CFSI 3.PLATF I 
FNUxCFSIA.PLATEI 
ALPH’CESI 5.PLATE ) 

BET.CFSIS.PLATEI 
GAM.CESIT, PLATE! 

RR«rFS(8*PLATE1 
r)«CFS(9. PLATE  I 
M»CFSI lO.PLATEI 

RSO«RR*RR 

ralrh*rr«alrh 

Rr>AM«RR*GAH 

f,MflrT»GAM*BFT 

CHIS0*« 0.9/BFT)»ALOGI ALPH/GAM) 

X2«-0.02»ALPH 
G?*h/ ( 1 .+FNU) 
r 1 1 •-?,*FNU) 

1 .-FNu) 
r‘^»ri»FNu 
r 

c SFT  UP  variablfs  for  Zone 
r 

10  S!GX*-S 

Sir,Y*OATPM  » 

S!G?«OATbI?> 

FFX.OATBt 

FFY*0ATB(41 

EEZ«nATB(5  » 

FPX.OATflJft) 

FPY-pATBI 71 
FPZ*OATB( 0) 

FTX. OATBJ9) 

FTY. OATBI lOI 
FTZsOATRdl) 
rB»oATB( 1 
XX*nATPI 1 3» 

R»OATP( 14  » 

RHsRR»B 

RBSO*RB»RB 

r 

r COMPUTE  strain  rates 

r 

nX»?*»lUN-UBN»/<0ELXj4.X-XR) 

0Y*o* 
r>z*o  • 

IFI  I ALPHA.FO.-nGO  TO  ZO 
OYs?.#<UN^UBN)/<  X4.Xfl>XN^XpN> 

IFI TALPHA.FO.l IPZ-OY 
?n  RHOf'OT»-pX-OY-nZ 
C 

c COMPUTE  Total  strain  increments  and  strains 

r 

OELTT«OFLT( 1 ) 

OFTX«nX»OELTT 

OFTY»f>Y«OFLTT 

DFTZ*OZ»nFLTT 

FTX*FTX>nFTx 

FTY«FTY4nFTY 

FTZ. FTZ^OFTZ 
C 

C SUPPOSE  INCREMENT  Is  ELASTIC  AnU  CHECK  FOR  YIELD 

C 
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r 

r 


NlNr»o 

Frx«FTX-FPK 

fFy.FTY-fpy 

FF2-FT2-FP2 

S!0X«C?»EFX4C3*<EEY+EEZ) 

Fji«S!r»x+siriY*siGZ 
FJ'^-FJl/1. 
bDX*SlGX-FJ3 
SDy«Sir»Y-FJ3 
SOZ«SIGZ.F J3 

F J?P«0*%»<  SOX«SOX^SDV*SOY^SOZ*SOZ I 
FJ?PR«S0RT<FJ2P) 

IFiFJl.GT.CRJGO  TO  30 
A»FJ1-CR 

FY«(RS0*FJ?P4A*A»/RBS0 

Fy»^ORT(FYI 

GO  TO  3? 

30  rF<FJUGT.CH!SOI  GO  TO  31 
FY«FJ?PR/(ALPH-GAm*EXP<BET«FJI ) ) 

GO  TO  3? 

31  FJ1*0» 

FY-lOO* 

CHECK  FOR  VlEunlNG 

3?  IF(FY.LE«1.  I GO  TO  60 
NrNC»(FY-l.>«?0» 

IF(nINC.GT.2000>  NINC8?000 
IF(NlNC,LT,n  NINC-1 


I 


C 

c 


r 

c 

c 


c 

c 

c 


RETURN  TO  ORIGINAL  STRAINS  AND  MAKE  INCREMENTS  SMALLER 

ETX«ETX-r>ETX 

ETY*ET Y-DETY 

FTZ«ETZ-nETZ 

F INC»NINC 

DETX-OFTX/FINC 

OETYsDFTY/FINC 

nETZ*f>FTZ/FINC 

USE  smaller  total  STRAIN  INCREMENTS 

00  66  INC«1  .NINC 
ETXsETX^DETX 

fty.fty^oety 

FTZ.FTZ+OETZ 

calculation  OF  YIELD  FUNCTION  VALUE 


‘kU  FFX*FTX-FPX 

FFY. ETY-FPY 

FFZ. FTZ-FPZ 

SIGX«C?»EEX^C3»*EEV>EEZ) 
SIGY.C?*EFY4C3*(EEZ>EEX) 
SIGZ«C?«EEZ4C3»<ECX4.EEY» 
F JI«5IGX^SIGY4SIGZ 
FJ3«FJl/3. 

SnX«SIGX-FJ3 

SnY«STGY-FJ3 

Sf>Z»SrGZ-FJ3 
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F J2P»0»6»I S0X*SDX4SDY»S0Y4S0Z«SDZ ) 

Fj?PR«S0RTIFj?P1 

IFIFJl.GT.CBIGO  TO  36 

A«rJT-CR 

F Y*|RS0*FJ?P4A»A|/RRS0 
GO  TO  38 


IFIFJ1*GT*CMIS0>  GO  TO  37 
FY«FJ?PR/(ALPH-GAM3EXPIBET«FJ1 » > 


GO  TO  30 

37  FJl-0. 
FY«100. 

38  CONTINUF 


r 


r 

ififv.le.1.0)  r,o  to  55 

r 

c calculate  plastic  strains 
c 

39  IF(FJl.GT.CB)  GO  TO  40 
YIELDING  ON  CAP  (F2) 

Al»l 

A2«?.»(FJl-CB)/ReS0 

BB3«.3.»A2*EXPI-D*I AX-XZ) )•( A242»*FT*BET»ll,-ALPH/bl )/ 
1 (ow*ii.+bet»(ralph-cb+xxi n 
IFLC-1 
GO  TO  so 

YIELDING  ON  FJ 


40  A1»FY/I2.»FJ2P) 

A2=2.»A1»FY»GMbET»FJ2PR«EXP(FJi*BET ) 

BB3*0. 

IFLG.O 

COMPUTE  gradients  WRT  stress 

So  GFX=41»SOX*A2 
GFY=41*Sr)Y+A2 
GFZ-A1»S0Z*A2 

CGFXiC2*GFX+C3*<GFY+GFZ> 

CGFYsC2»GFY+C3»(GFZ»GFX) 

CGFZ«CZ*GFZ+C3«IGFX+GFY) 

BBI=CGFX*dETX+C&FY*OETY+CGFZ»DETZ 

BB2=CGFX*GFX*CGFY«GFy+CGFZ«GFZ 

nLAM=BBl/(eB2-BB3l 

COMPUTE  new  plastic  STRAINS 

DEPX»OLAM*GFX 

depy*dlam»gfy 

OEPZ«OLAM*GFZ 

FPX=EPX*OEPX 

fpy*epy+oepy 

EPZ.EPZ+OEPZ 

EP1=EPX+EPY+EPZ 

NINNERoO 

ARG.l.O+EPl/W 
IFIARG.LE.0»0)  go  To  112 

xx^xz+ialociwo-aepi/wi  i/o 

XRA.XXaRALPH 

CALL  COMPC(CB*XRA.RgaM.B£T I 

P»lrfl-XX)/RR 

RB*rB-XX 

RBSQ<RB*RB 

GO  TO  113 

112  EEX.FTX-FPX 
EEY.ETY-EPY 

eez.etz-epz 

SIGX«C2»EFX+C3» I EE Y+EEZ ) 

SIGY.C2»EEY+C3*IEEZ*EEX) 

SIGZ*C2»E€Z*C3*IEEX+EEYI 

FJl.SIGX+SIGY+SIGZ 

CB«FJ1 

B>AlPH-GAM*EXP(BET«FJ1 I 

RB«RR»B 

RBSo«RB>RB 

XX»CB-R8 

GO  TO  1 14 

calculate  yield  function  with  new  values  OF  plastic  strain 

113  FFX.ETX-EPX 

eey.ety-epy 
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r\  n rt  r\  rx'^ 


r 


♦C3*«Eev+Et2» 

Sir>Y*C?*EEY^C3»lEEZ>EEX» 

SIGZ«C?*EEZtC3»(tEX*EEYI 
FJl*SlGX^SirfY^SlG2 
1 14  FJ3*FJl/3, 

SnX«S|GX-FJ3 
snY.sinY-.Fj3 
SnZ»S!GZ-F J3 

F J?P«0*^*t  S0X*S0X4Sf>Y»S0Y4SD2«S0Z  » 

FJ?PR«S0RTtFJ2P) 

IF*FJl.r>T.CB)  GO  TO  46 

r 

A*FJ1-C8 

FY«(  RSO*FJ2P«*-A»Ai/RBSa 
GO  TO  41 

48  IF(FJ1.GT.CHIS0»  GO  TO  49 

FY«FJ>PR/IALPH-GAM*EXP(PET#FJI ) ) 

GO  TO  41 

40  rji*n, 

FY«100, 

41  rONTlNUF 

C 

C 0.99.Lt.FY.LE.  1.01  IMPLIED  SATISFACTORY  STRESS  PoINl. 

IFI*-Y,GT.1.01  »G0  TO  39 
IFIFY.GT.0.991GU  TO  SS 

THIS  PORTION  OF  THt  ITERATION  PREVENTS  OSCILLATIONS  APOUT  A POINT 
ON  the  yield  surface,  plastic  deformation  has  OCCURRED  BuT 
FY  .lT.0.99.  reduce  PLASTIC  STRAINS  TO  APPROACH  YIELD  SURFACE 
FROM  EXTERIOR. 

NINNFR»NINNFR+1 
IF1NINNER.GT.S0)  GO  TO  5S 
FACT*-0. 1 
FPX*EPX^FACT*0EPX 

EPY.EPY-*-FACT*DEPY 
EPZ»EPZ^FACT»nEPZ 
DFPX«0.9»DEPX 
DEPV*0.9*r>FPY 
nFPz«n.9*0EPz 
FP1«FPX+FPY*FPZ 

ARG«1 .+FP1/W 
IFIARG.LE.0.01  go  to  112 
XX«XZ>1 ALOGI 1 .O^EPl/W) J/D 

xra.xx^ralph 

CALL  COMPClCBfXRA.RGAM.BET) 

B*<CR-XX)/RR 
RBsRRvR 
RRSO*RR»Rfl 

GO  TO  in 

AS  rONTINlIF 
AO  rONTINUF 

( STORE  NFW  variables 

c 

SN«-SIGX 

PN*.SlCiX.SIGY-SlGZ 

ZN*SIGX-SIGY 

FN-P 

rN«ro 

OATo I 1 ) «SIGY 

0ATR1?)»SIGZ 

r>ATnl3)aFFX 

r)ATAU)*EFY 

nATo(S»*FFZ 

nATR(A)«FPX 

nATR(7)aEPY 

DATR(8)«EPZ 

nATRl9)«FTX 

datrmo)*ety 

OATRllDaFTZ 

OATRini-CB 

OATR 1 1 3 ) aXX 

OATo 1141*0 

nATomiaFY 

OATolm^FJ! 

OATpI 17)*FJ2PR 
OATrI 18)»FP1 
PFTDRN 
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t 


I 

! 

f 

( 


FNTRV  STIN? 

IFiMfwPLATlSrtfttSlO 
SftO  •LPH«rFSt S*PLATF » 

PFT*rF<i(6*PLAtE  » 
r,AH«CFS(  7*PLATFi 

QR*rF<( RtPLATF)  \ 

r>«fFSlQ*PLATFl  ^ 

RAtPH«RR*ALPM 
RGAm«RR«GAM 
1 

XK*-0.0?»ALPH 
XHA*XX>PALPM 

call  compckb*xra*roam*bet > 

^ATpC 

r>ATR(nJ«XX 
R« I rR-XX ) /RR 

OATr I 1A ) *r 
rfturn 

f JO  oaTqi i?j»rR 
nATRll^)»XX 
0AT«(  I41«p 
RFTuRN 
ENn 


subroutine  C0MPC(C*C1*C2*C3» 
r 

C NFWTnNS  MFTHOO 

r 

COLn»C 

NC*0 

10  At«r?»FxPic3*r ) 

A?«ri-Al-C 

Al»r3«Al4l« 

r*r4A2/A3 

IF(ABSCC3«IC0LD-CM  .lT.O.OOOI  IGO  TO  30 
NC«NC^1 

|FI(«4C.GT*20>GO  to  20 

roLn.f 

GO  TO  10 

70  WRITFI6M0001C 

1000  FORmaTIIH  ♦AX*16HC  NOT  CONVERGING  •E10*3) 
STOP 

*^0  ific.gt.o#ok*o* 

RETl/RN 

END 


J 

1 
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ABBREVIATIONS  AND  SYMBOLS 


C 

[C] 

^ijkt 

D 

E 

K 

P 

R 

W 

X 

Z 

b 

Kj 

dA 

e p 

{e®} 

t 

a,  3,  Y 

P 


e 
w 

V 

P 

{a} 


e p 
M , 


Jj  value  at  center  of  ellipse  of  cap  model 

constitutive  matrix 

compliance  components 

cap  model  parameter 

Young's  modulus 

first  invariant  of  stress  tensor  = o - - 
second  invariant  of  stress  deviator  tensor  = ^ 
bulk  modulus 
pressure  = -0^/3 

ratio  of  major  to  minor  axis  of  ellipse  of  cap  model 
cap  model  parameter 

value  at  intersection  of  cap  and  J^-axis 
J value  at  initial  position  of  cap 
VST  value  at  center  of  ellipse  of  cap  model 
total  strain  rate  or  increment 

elastic  strain  rate  or  increment 

plastic  strain  rate  or  increment 
scalar  function 

total,  elastic,  and  plastic  components  of  strain  tensor 

elastic  strain  vector 
yield  functions 
time 

cap  model  parameters  which  define  ideal  yield  surface 
Kronecker  delta 

volumetric  plastic  strain 

total,  elastic,  and  plastic  volumetric  compression 
Poisson's  ratio 
mass  density 
stress  components 

deviator  stress  components 

stress  vector 

plastic  potential  function 


J 
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